TWO-GRID METHODS FOR SEMILINEAR INTERFACE PROBLEMS 

MICHAEL HOLST, RYAN SZYPOWSKI, AND YUNRONG ZHU 

Abstract. In this article we consider two-grid finite element methods for solving 
semilinear interface problems in d space dimensions, for d = 2 or d = 3. We first 
describe in some detail the target problem class with discontinuous diffusion coeffi- 
cients, which includes problems containing sub-critical, critical, and supercritical non- 
linearities. We then establish basic quasi-optimal a priori error estimate for Galerkin 
approximations. In the critical and subcritical cases, we follow our recent approach to 
C^) controling the nonlinearity using only pointwise control of the continuous solution and 

CN a local Lipschitz property, rather than through pointwise control of the discrete solution; 

this eliminates the requirement that the discrete solution satisfy a discrete form of the 
^^ maximum principle, hence eliminating the need for restrictive angle conditions in the 

^ j underlying mesh. The supercritical case continues to require such mesh conditions in 

. order to control the nonlinearity. We then design a two-grid algorithm consisting of a 

coarse grid solver for the original nonlinear problem, and a fine grid solver for a lin- 
earized problem. We analyze the quality of approximations generated by the algorithm, 
^h and show that the coarse grid may be taken to have much larger elements than the fine 

7 j grid, and yet one can still obtain approximation quality that is asymptotically as good 

as solving the original nonlinear problem on the fine mesh. The algorithm we describe, 
and its analysis in this article, combines four sets of tools: the work of Xu and Zhou 
on two-grid algorithms for semilinear problems; the recent results for linear interface 
problems due to Li, Melenk, Wohlmuth, and Zou; and recent work on a priori estimates 
for semilinear problems. 
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1. Introduction 

In this article, we consider a two-grid finite element method for semilinear interface 
problems with discontinuous diffusion coefficients. One of the primary motivations of 
this work is to develop more efficient numerical methods for the nonlinear Poisson- 
Boltzmann equation, which has important applications in biochemistry and biophysics. 
However, the theory and techniques are applicable to a large class of semilinear inter- 
face problems, including problems with critical (and subcritical) nonlinearity arising in 
geometric analysis and mathematical general relativity. 

In order to achieve our goal of exploiting two-grid-type discretizations, our first task 
is to more completely develop a basic quasi-optimal a priori error analysis for Galerkin 
approximations of semilinear interface problems. The main challenge comes from the 
loss of global regularity for interface problems (cf. [3, 17]). There has been much work 
on finite element approximation of the linear elliptic interface problem. For example, 
in [3] an equivalent minimization problem was introduced to handle the jump interface 
condition; this problem was then solved using finite element methods. Subsequently, fi- 
nite element approximation of the elliptic interface problems was analyzed using penalty 
methods in [7], and optimal rates in the H 1 and L 2 norms were obtained by appropri- 
ately choosing the penalty parameter. Optimal a priori error estimates for linear interface 
problems in the energy norm (i.e., a weighted H 1 norm) is given in [25]. In [31], sub- 
optimal error estimates of order 0{h\ \ogh\ 1 / 2 ) in the H 1 norm was obtained for 2D 
linear interface problems using standard finite element techniques. Similarly, in [11] it 
was shown that for C 2 interfaces in 2D convex polygonal domains T, the linear FEM 
approximation u h has suboptimal standard error estimates of orders 0(h\ log h\ 1 / 2 ) and 
0(h 2 \ log/i] 1 / 2 ) in H 1 and L 2 norms respectively. By using isoparametric elements to 
fit the smooth interface, optimal error estimates were obtained in [27] for 2D interface 
problems. These results have been generalized to higher-order finite elements approxi- 
mation in [21]. There are also other approaches for dealing with linear elliptic interface 
problem; for example, immersed interface finite element methods based on Cartesian 
grids (cf. [22]), mortar finite element (cf. [20]), and Lagrange multiplier methods using 
non-matching meshes (cf. [13]). 

Less work has been done for nonlinear interface problems. For smooth coefficients un- 
der quite strong (global) regularity assumptions, quasi-optimal error estimates were ob- 
tained by [32, 34]. Due to the loss of global regularity for interface problems (cf. [3, 17], 
see also [23, 24] for regularity of linear interface problems), these analysis techniques are 
not applicable here. Recently, Sinha and Deka [28] studied linear finite element approx- 
imation of semilinear elliptic interface problems in two dimensional convex polygonal 
domains. Under assumptions that the mesh resolves the interface, and that the nonlinear 
function &(£) satisfies 

\V(Z)\<C\Z\, and|6"(0I^C V£ e R, 

they showed optimal error estimates in the H 1 norm using the framework of [8], together 
with the results from [11]. 

In this paper, we use a more natural approach for semilinear interface problems which 
can be applied to a somewhat different but larger class of nonlinear problems than [28]. 
For ease of exposition, we assume that the triangulation resolves the interface, although 
this assumption may be weakened. The first step is to derive both continuous and dis- 
crete a priori L°° bounds for the continuous and discrete solutions in order to control 
the nonlinearity. While continuous L°° bounds are fairly standard under quite general 
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assumptions on the nonlinearity (cf. Assumption 2.2), discrete a priori L°° bounds re- 
quire additional mesh conditions on the triangulation (cf. Assumption 3.1). Based on 
a priori L°° control of the continuous and discrete solutions, we derive optimal a pri- 
ori error estimates in both the H 1 and L 2 norms, with the help of a Local Monotonicity 
assumption on the nonlinearity. A similar approach was used in [9, 14] for the Poisson- 
Boltzmann equation. We note the mesh conditions play a key role in obtaining discrete 
maximum/minimum principles (cf. [18, 15, 16, 30]). However, when the nonlinearity 
satisfies subcritical or critical growth conditions, and has some type of monotonicity, we 
have been able to derive quasi-optimal a priori error estimates directly, without using 
discrete maximum principles, and hence without the need for any mesh angle condition 
assumptions [4]. 

Finite element approximation of semilinear interface problems results in the need to 
solve system of nonlinear algebraic equations, and the number of unknowns in these 
systems can be extraordinary large in the case of three or more spatial dimensions. The 
most robust and efficient approach for solving these types of nonlinear algebraic systems 
has been repeatedly shown to be some variation of damped inexact Newton iteration, 
which consists of an inner-outer iteration: an inner loop involving repeated linear solves, 
together with any outer loop involving a damped/inexact correction step. See for exam- 
ple [5, 6, 26], and also [1] for an application to nonlinear interface problems. The basic 
approach involves the solution of a linear system on the fine mesh at each Newton step. 
However, the two-grid algorithm proposed in [2, 32] takes another approach, which con- 
sists of a coarse grid solver for the original nonlinear problem, and a fine grid solver only 
involving for a linearized problem, which is effectively a one-step Newton update of the 
solution. The benefit of using this two-grid idea is that it significantly reduces the overall 
computation cost, since we only need to solve the nonlinear problem on a coarse grid, and 
we can solve the linear problem on the fine grid by using standard multigrid/multilevel 
methods for optimal complexity. The central question concerning the two-grid method 
is to how choose the coarse grid problem; in other words, how coarse can one make the 
nonlinear problem discretization, but still achieve nearly optimal approximation proper- 
ties if solving the full nonlinear problem on the fine grid. Based on a priori H 1 and L 2 
error estimates for semilinear interface problems, in this paper we show that the basic 
framework developed in [32, 33] allows us to establish, both theoretically and numeri- 
cally, that one may choose a coarse grid with much larger mesh size than the fine grid in 
the case of semilinear interface problems. 

The main contributions of this paper are as follows: 

(1) We give a complete finite element error analysis for semilinear elliptic interface 
problems, under weak assumptions on the nonlinearity; this includes establishing 
quasi-optimal a priori energy, L 2 and L 4 error estimates of the finite element 
approximation. 

(2) We also provide a practical approach to efficiently solve the resulting nonlinear 
algebraic problem by two-grid algorithms, reducing the solution of the original 
nonlinear system of equations on the fine grid to the solution of a nonlinear prob- 
lem on a coarse grid having much fewer degrees of freedom, together with the 
solution of a linear problem on the fine mesh. We note that the resulting linear 
interface problem can be efficiently solved by PCG algorithms using multilevel 
or domain decomposition preconditioners (cf. [35, 36]). 

The remainder of the article is organized as follows. In Section 2, we introduce the 
basic notation and the model problem. We also establish continuous L°° bounds for the 
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solution under very weak assumptions on the data and the nonlinearity. In Section 3, 
we establish quasi-optimal error estimates for the finite element approximation, by first 
deriving discrete a priori L°° bounds. In Section 4, we describe the two-grid algorithm, 
and give an analysis of the approximation properties of the algorithm. In Section 5, we 
give some numerical experiments to support our theoretical conclusions. 

2. Semilinear Interface Problems 

Let ft C R d be a Lipschitz domain with d > 2, with an internal interface T dividing it 
into two open disjoint subdomains fti and ft 2 , so that ft = fti U T U ft 2 . For ease of expo- 
sition, we assume ft x and ft 2 are two non-overlapping polyhedral/polygonal subdomains. 
We then focus on the following semilinear elliptic equation: 

- V • (DVu) + 6(x, u) = in ft, u\ dQ = 0, (2.1) 

with the jump conditions on T: 

du 



[u] =0, and 



on 



on T, (2.2) 



where [u] := ui\ r — u 2 \ r and [D^] := i^if^- + D 2 ^ with n^ representing the unit 
outer normal on ft^. Here ui (i — 1, 2) stands for the restriction of u on ft^. We assume 
that the coefficient D = D{x) : ft — >• ]R (ixci is symmetric and piecewise constant on each 
subdomain, i.e., D = Di in ft^, and D G L°°(ft) satisfies 

m\i\ 2 ^ i T D(x)i < M|£| 2 , V£ GR d , x G ft, (2.3) 

for constants m, M > 0. 

In working with the solution and approximation theory for (2.1)-(2.2), we will employ 
standard notation for the function spaces, norms, and other objects that will be relevant. 
For example, given any subset G C R d , we denote as L P (G) the Lebesgue spaces for 
1 ^ p ^ oo, with norm || • ||o, p ,g- We denote the Sobolev norms as ||v|| 5jPj g — IM|w s >p(G) 
for the spaces W S > P (G), with W S > 2 (G) = H S (G) when p = 2. For any two functions 
v G L P (G) and w G L q (G) with p, q > 1 and 1/p+l/q = 1, we denote the pairing 
(v,w)g '•— J G vwdx. For simplicity, when G = ft, we omit it in the norms/pairings. 
We will also denote as H s (Qi U ft 2 ) the space of functions u such that u\q. G H s (Qi) 
for i = 1, 2 and 5 > 1, endowed with the norm 

Hll2 ||||2 ||||2 

ll^llif s (niun 2 ) := IMI# s (^i) + \\ u \\h s (q 2 )- 

We will use the notation x\ < y\, and x 2 > y 2 > whenever there exist constants Ci, C 2 
independent of the mesh size h and the coefficient D or other parameters that xi, x 2 , 
yi and y 2 may depend on, and such that x\ ^ Ciyi and x 2 > C 2 ?/ 2 . We also denote 
x ~ y as Cix ^ y ^ C 2 x. Without confusion, we will also write &(£) := &(#,£) and 
&'(£) := 96(x, £)/^£ f° r simplicity. 



Remark 2.1. Afote ^/itftf more general interface conditions D^- = grfor some given 

function gr G L°°(r) and non-homogeneous Dirichlet data u\oq = g can be easily 
treated using our results here, due to the observation that one may split the equation 
into two sub-problems. The first sub-problem is a linear elliptic interface problem, and 
the second sub-problem is a nonlinear elliptic problem (2.1) with homogeneous Dirichlet 
boundary condition. More precisely, let u = u l + u n , where u l G Hg(Q) satisfies the 
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linear elliptic interface problem: 

-V • (DVu 1 ) = in CI 

u l \dn = 9: and 

while the nonlinear part u n is the solution to the (homogeneous) semilinear equation 

- V • (DVu n ) + b(u n + v!) = in ft, 

with the interface condition (2.2). On the other hand, the treatment for the linear in- 
terface problem (2.4) is standard; cf [11, 21]. Therefore, without loss of generality we 
focus on (2.1) with homogeneous interface conditions (2.2). 

The weak form of (2.1) reads: Find u e Hq(Q) such that 

a(u, v) + (b(u),v) = 0, W e H^Q), (2.5) 

where a(?i, v) := J n L>Vt£ • Vvdx. By the assumption (2.3) on the coefficient D, the 
bilinear form a(u,v) in (2.5) is coercive and continuous, namely, 

ra||V^||o 5 2 ^ a(u,u), a(u,v) ^ M||V^||o,2||V^||o,2 5 \/u,v e i?o(fl), (2.6) 

where < m ^ M < oc are constants depending only on the maximal and minimal 
eigenvalues on D and the domain Q. The properties (2.6) imply the semi-norm on Hq (Q) 
is equivalent to the energy norm ||| • ||| : Hq(Q) —> M, 

||k||| 2 = a(u,u), mllV^Hno ^ IIMII 2 ^ MllV^Hno- (2.7) 

MINI \ ' / ' II IIU.^illllll II I I \J * £j x * 

A priori L°° bounds for any solution to the continuous problem play a crucial role in 
controlling the nonlinearity. The following weak assumption on the nonlinearity allows 
for a large class of nonlinear problems containing both monotone and non-monotone 
nonlinearity: 

Assumption 2.2. b : Q x R —> R is a Caratheodory function, which satisfies the barrier- 
sign conditions in its second argument: there exist constants ft,^G| with a ^ f3, such 
that 

b(x,£)^0, Vf^/3, a.e. in Q 
b{x ) £) ^0, V<^ ^ a ) a.e. in Q. 

We have the following theorem based on the Assumptions 2.2: 

Theorem 2.3 (A Priori L°° Bounds). Let the Assumption 2.2 hold. Let u e H^{Vt) be 
any weak solution to (2.5). Then 

u^u^u, a.e. in fi, (2.8) 

for the constants u and u defined by 



u := max < /?, sup g(x) > , u \— min < a, inf g(x) > , (2.9) 

[ xedn J [ ^<^ J 

where a ^ /3 are £/ie constants in Assumption 2.2. 
Proof. To prove the upper bound, let us introduce 

(j) = (u — u) + = IRdiX{u — u, 0}. 

By the definition of % it follows (cf. [29, Theorem 10.3.8]) that (j) e #o(^) and > 
a.e. in fi. Taking i> = in (2.5), we have 

a(u, (p) = a(u -%</>) = a((f>, (j)) ^ m\\ V0||^ 2 . 
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This implies that 

m||V0||g j2 <a(Tx,0) = (-6(^),0)<O, 

since — b(u) ^ a.e. in the support of cj). Hence, || V0||o,2 = which yields = 0. 
Therefore, the upper bound of (2.8) holds. 
Similarly, we introduce 

ij) = (u — u)~ = min{u — u, 0}. 

It is obvious that i/j e Hq(Q) can be used as a test function in (2.5). Moreover, ip ^ 
a.e. in Q, and Assumption 2.2 implies —b(u) > on the support of i/j. Therefore, 

which implies r/j = as before. This proves the lower bound of (2.8). □ 

To conclude this section, we give the nonlinear Poisson-Boltzmann equation as an 
example, which is one of our main motivation for this work. This equation has been 
widely used in biochemistry, biophysics and in semiconductor modeling for describing 
the electrostatic interactions of charged bodies in dielectric media. 

Example 2.4. The regularized Poisson-Boltzmann equation reads: 

- V • (sVu) + k 2 sinh(^) = 0, in Q 

[u]r = 0and[e^] r = <?r, onY (2.10) 

u\dn = 9, ondVt, 

where gr € L°°(r) is a function defined on Y arising from regularization of pointwise 
charges in the molecular region (see [9, 14] for detailed derivations). Here the diffu- 
sion coefficient e is piecewise positive constant e|^ 1 = e\ and s\q 2 = s 2 , where Qi is the 
molecular region, and Q 2 is the solution region(see Figure 2.4 for example). The modified 




Debye-Huckel parameter k 2 is also piecewise constant /^ 2 (x)|^ 1 = and n 2 {x)\^ 2 > 0. 
The Dirichlet condition u\qq = g are imposed on the boundary dQ. We note that equa- 
tion (2.10) can be reduced to (2.1) by splitting it into linear and nonlinear components 
as described in Remark 2.1, see [\A]for more details. Obviously, the Assumption (Al) is 
satisfied for (2.10). 

3. Finite Element Error Estimates 

We now discuss some error estimates on the finite element discretization of (2.5) which 
will play a key role in the two-grid analysis. Given a quasi-uniform triangulation % of 
H, we denote by Vh(g) C Hg(Q) the standard piecewise linear finite element space 
satisfying the Dirichlet boundary condition. For simplicity, we denote V h := V/^O). For 
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ease of exposition, we assume the triangulation Th resolves the interface Y. Then finite 
element approximation of the target problem (2.1) reads: Find u h G 14 such that 

a(u h) v) + (b(u h ),v) = 0, Vve V h . (3.1) 

The following theorem shows that under appropriate mesh condition on Th, the discrete 
solution u h of (3.1) satisfies a priori L°° bounds (as does the continuous solution u due 
to Theorem 2.3). More precisely, assume the triangulation Th satisfies 

Assumption 3.1. Let fa and fa are the basis functions corresponding to the vertices X{ 
and Xj, respectively. We assume that 

X 0, Vi ^ j, (3-2) 



Under this assumption, we can obtain the following a priori L°° bound of the discrete 
solution u^ 

Theorem 3.2. Let b satisfy the Assumption 2.2, and Th satisfy the Assumption 3.1. Then 
the discrete solution Uh € Vh(g) to (3.1) satisfies 

u^Uh^u, a.e. in fi, (3.3) 

where u and u are defined in (2.9). 

The idea of the proof of (3.3) is the same as in Theorem 2.3. However, in the discrete 
setting, for a given r G R the truncated functions (u h — r) ± are usually not in Vh. 
Thus, they can not be used as test functions in (3.1). Instead, one can employ the nodal 
interpolation of these functions as test functions. In particular, given any constant r, we 
denote 

N 

K - r ]^ : = ^2( u h(xi) - r^fa, 

2 = 1 

where TV is the total number of degree of freedoms, and fa (i — 1, • • • , N) is the nodal 
basis function at the vertex xi. While this does produce proper test functions, it unfortu- 
nately introduces mesh conditions such as Assumption 3.1 into the analysis. 

Proof of Theorem 3.2. To prove the upper bound of (3.3), define a test function (p^(x) := 
[uh(x) — u] + . It is obvious that fa^ e Vh and supp(0^) is the union of of the macro 
elements for the vertices such that u h {xi) > u. Therefore, it satisfies 

aK,^) + (5K),^) = 0. (3.4) 

For the diffusion term in (3.4), we notice that 



a ( u h,(/>h) = a (0h^h)+ a (i U h 



U 



,J h. 



a m > ■ 



^h) + y]( u h(xj) ~ u ) ( u h(xj) ~ u) + a(fa, fa) 

> a(<^+)^m||V<^|lo, 2 , 

where we used the Assumption 3.1 in the third step, and used (2.3) in the last step. Thus 
we obtain that 

H|V<^|lo,2 < <u h At) = {-b{u h )At) < 0. 

This implies that fa^ = 0, and hence u h ^u a.e. in Q. The proof of the lower bound is 
similar, and so we omit the detail here. □ 
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Remark 3.3. Note that Assumption 3.1 requires certain angle condition on the triangu- 
lation Th- This condition is crucial in proving the discrete maximal/minimal principle (cf. 
[12, 18, 15, 30]). However, in case that b satisfies critical/ subcritical growth condition, 
namely, there exists some constant K > such that 

\b {n \0\< K , V£eM (3.5) 

where n is an integer satisfying n < oc when d = 2 and n ^ (d + 2)/{d — 2) when 
d ^ 3, we are able to show the quasi-optimal error estimate directly, without using 
Assumption 3.1; see [A] for more detail. 

A priori L°° bounds (Theorem 2.3 and Theorem 3.2) play crucial roles in controlling 
the nonlinearity, ensuring that the nonlinearity b has a certain "local Lipschitz" property. 
This property in turn is used to establish quasi-optimal error estimates for the finite ele- 
ment approximations. For this purpose, let us make the following additional assumption 
on b : 

Assumption 3.4. b is locally monotone, namely, 

b'(0>0, ^e[u,u], (3.6) 

where u, u are the barriers defined in (2.9). 

Without loss of generality, in the remainder of the paper, we let the Dirichlet data 
g = 0. With the help of the a priori L°° bounds of u and u h in Theorem 2.3 and Theo- 
rem 3.2 respectively, we are able to establish the following quasi-optimal error estimate. 

Theorem 3.5. Let b satisfies the Assumptions 2.2 and 3.4, and Th satisfy the Assump- 
tion 3.1. Let u G Hq (Q) n H s (Qi U Q2) f or some s > 1 and Uh G 14 be the solution 
to (2.5) and the discrete solution to (3.1), respectively. Then the following quasi-optimal 
error estimate holds: 

\\\ u ~ u h\\\ < inf \\\u - v\\\ < ^"^Mlij^ftiUfia). (3.7) 

vev h 

Proof. By Assumptions 2.2 and 3.1, Theorem 2.3 and 3.2 give a priori L°° bounds on u 
andt^: 

U ^ U ) Uh ^ u 

for the constants u^u defined in (2.9). This implies that 

b(u) - b(u h ) = &'(0(u - u h ) ^C L \u- u h \, (3.8) 

where C L = sup^^j ||&'(0 llo,oo * s a constant depending only on &, u and u. Subtracting 
equation (3.1) from (2.5), we have 

a(u - u h , v) + (b(u) - b(u h ),v) = 0, \/v G V h . 

By using this identity, we obtain 

\\\u — Uh\\\ 2 = a(u — Uh, u — Uh) 

= a(u — Uh, u — v) + a(u — u^ v — Uh), Vv G 14 

= a(u -u h ,u-v) + (b(u) - b(u h ),u h - v) 

= a(u -u h ,u-v) + (b(u) - b(u h ),u - v) - (b'(g)(u - u h ),u - u h ) 

^ l U - U h\\u - V\ + C L \\U - U h \\o^\\u - ^||o,2, 
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where in the last inequality, we used Cauchy-Schwarz inequality, the Lipschitz property 
(3.8) of b and the Local Monotonicity (3.6) from Assumption 3.4. Then by Poincare 
inequality we have 

WW ~ U h\\\ 2 <: \\U - U h \\\\\\u - V\\\ + C L Cp\\V{u - ^)|| ,2||V(^ - V)|| 0j 2 

where C P is the Poincare constant. Thus we obtain 

l u ~ u hl ^ Ilk - v \t 

Therefore, we have proved the first inequality in (3.7), since v G Vh is arbitrary. The 
second inequality in (3.7) follows by standard interpolation error estimates; cf. [21, The- 
orem 3.5]. □ 

To conclude this section, let us try to derive L 2 error estimates: \\u — Uh\\o,2> by using 
duality arguments. To begin, introduce the following linear adjoint problem: Find w G 
i?o(fi) such that 

a(v, w) + (b'(u)v, w) = (u- u h , v), v G Hq(Q). (3.9) 

We assume that the linear interface problem (3.9) has the regularity 

IMItf T (ftiun 2 ) ^ C \\ u ~ ^llo,2 (3.10) 

for some r > 1. The regularity assumption (3.10) is also called "r-regularity" in [21, 
Assumption 4.3], which is quite natural for linear interface problems. Along with (3.9), 
let us also introduce the finite element approximation w h G Vh which satisfies: 

a(v h ,w h ) + (b'(u)v h ,w h ) = (u-u h ,v h ), v h G V h . (3.11) 

Then by standard finite element approximation theory for the linear interface problem (3.9) 
(cf. [11, 21]), we have the following error estimate: 

III™ - w h\\\ < ^"iHli^Cftiuna) ~ ^ r_1 ||^ - ^||o,2, (3.12) 

where in the second inequality we have used the regularity assumption (3.10). We then 
have the following L 2 error estimate for u h : 

Theorem 3.6 (L 2 Error Estimate). Let b satisfy Assumptions 2.2 and 3.4, and let % 
satisfy Assumption 3.1. Let u G Hq(Q) PI H s (Qi U Q2) with s > 1 be the solution 
to (2.5), and let u^ he the solution to (3.1). Suppose that the dual problem (3.9) satisfies 
the r-regularity (3 .10) for some r > 1. Then 

\\u - u h \\ ,2 < C(h T ~ l + h-^WViu - u h )\\ 0t2 , (3.13) 

where C is independent ofh. 

Proof. Without loss of generality we assume b"{x) ^ 0. By taking v = u — u h in (3.9) 
we obtain that 

11^ - ^llo,2 = a ( U ~ U h, W) + (b\u)(u - U h ),w) 

= a(u -u h) w - w h ) + (b'(u)(u - u h ),w - w h ) 

+ a(u - u h) w h ) + (b'(u)(u - u h ),w h ). (3.14) 

To bound the second term in (3.14), we use the L°° bound of u (cf. Theorem 2.3) to 
obtain 

(b'(u)(u - U h ),W - W h ) ^ ||&'(tz)|| ,oo||u - ^||o,2|k - W/i||o,2 

< \\\7(u - ^)|| ,2||V(^ - W h )\\o,2, 
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where we have used Poincare inequality for (u — u h ) and (w — w h ) in the last step. To 
deal with the last two terms in (3.14), notice that u h G 14 is the solution to the discrete 
semilinear problem (3.1), we have 

a(u - u h) w h ) = ~(b(u) - b(u h ),w h ). 

Thus by Taylor expansion, we have 

a(u - u h: w h ) + (b'(u)(u - u h ),w h ) = -(b(u) - b(u h ) - (b'(u)(u - u h ),w h ) 

= -{b f \x)(u-u h ) 2 ,w h ), 

where x satisfies that u ^ x( x ) ^ ^ a - e - i n ^ due to the a priori L°° bounds of u and u h 
by Theorem 2.3 and Theorem 3.2. Therefore, by Holder inequality the last two terms in 
(3.14) can be bounded as: 

a{u-u h) w h ) + (b'(u)(u-u h ),w h ) ^ - ||&''(x)||o,oo||u - u h\\l, p *\\ w h\W q * 

< ||V(TX-^)||o, 2 |k/i||o, 9 * 3 

where we choose p* = 6 for d = 3 and p* > 4 when d = 2, and 4" + -^ = 1. In the 
last inequality, we have used the Sobolev embedding \\u — Uh\\o, P * ^ ||V(m — iOi)||o,2- 
Therefore, we obtain 

11^ - U h\\l,2 ^ (ll V (^ - W/0Ho,2 + ||V(iX - ^)||o,2||^||o,g*)ll V (^ ~ ^)l|o,2- (3.15) 

Now in (3.15), let t^ e 14 be the solution to (3.11). Then the estimate for the quan- 
tity \\V(w — w h )\\ ^ is readily available from (3.12). It then remains to estimate the 
term \\V(u — ^/i)||o,2||^/i||o,g*- To estimate ||^||o,g*» notice that by the choice of p* 9 
1 < g* < 2. Then by Poincare inequality and coercivity of a(-, •) 

IWkg* ^ \\ W h\\o,2 < \\Vw h \\ 0j 2 ^ Ilk/ill- 

By Assumption 3.4 on 6, we have 

\\\w h \\\ 2 = a{w h) w h ) ^ a{w h) w h ) + (b'(u)w h ,w h ) 

= (u-U h ,W h ) ^ ||w-^||o,2||^||o,2 



< 



\v> ~ U h \\o^\\\W h \\ 



Thus, this inequality implies that 
Therefore, we obtain 



Wh\\o,q* ^ llk/zll ^ \\U-U h \\ ,2. 



\\V(u - u h )\\o^\\w h \\ ^ < h s 1 \\u\\ H s(n 1 \jn 2 )\\u - u h \\ ^- (3.16) 

Combining inequalities (3.15), (3.12) and (3.16), the inequality (3.13) then follows. This 
concludes the proof. □ 

Remark 3.7. Theorem 3.6 implies, in particular, that 

\\u - u h \\^ 2 % /&*||V(tx - ^)||o,2, (3.17) 

where t = min{5,r} — 1. Combining with the quasi-optimal error estimate (3.7), we 
obtain the following L 2 error estimate: 

\\u - u h \\o, 2 < ^+( s - 1 )||^|| iJS(0iUa2) . (3.18) 

Remark 3.8. Similar to Remark 3.3, ifb(-) satisfies the critical/ subcritical growth con- 
dition (3.5), then the conclusion of Theorem 3.6 holds without the Assumption 3.1. We 
refer to [4] for the details. 
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4. Two-grid algorithms 

We now consider a two-grid algorithm (cf. [32]) to solve the finite element discretiza- 
tion (3.1) numerically. Let Th be a quasi-uniform triangulation with mesh size H, and 
Th with mesh size h < H is a uniform refinement of Th • We assume the triangulations 
satisfy Assumption (3.1). The algorithm consists of an exact coarse solver on a coarse 
grid Th, and a Newton update on the fine grid Th- In what follows, we will denote u h) u H 
as the exact finite element solutions to (3.1) on the grids Th and Th, respectively. For 
simplicity, let us denote 

(F(u),x) '•= ft (^x) + (&H>x)> 
and its linearization 

(F'{u)v, x) '= Q>(v, x) + 0>'(u)v, x)- 
The two grid algorithm considered in this paper is as follows: The Algorithm 1 solves 

Algorithm 1: u h = TwoGrid (T H , %) 

1 Find uh G V h such that 

(F(u H ),v H ) = 0, Vv H eV H l 

2 Find u h G Vh such that 

(F\u H )u h , v h ) = {F'(uh)u h , v h ) - (F(u H ),v h ) Vv h G V h ] 

the original nonlinear problem on the coarse grid 7#, and then performs one Newton 
iteration on the fine grid. 

Fix any Xh G V h , let rj(t) := (F(u H + t(u h - u H )),Xh)- Then by Taylor expansion, 
we have 

= (F(u h ),Xh)=v(l)=v(0)+v\0)+ I v"(t)(l-t)dt 

Jo 

= (F(u H ), X h) + (F\u H )(u h -u H ),Xh) + [ v"(t)(l-t)dt. 

Jo 

Notice that by direct calculation, the remainder term, denoted by R{u H) u h) Xh), has the 
following form: 

R(u H ,u h) Xh) '•= / r)"(t)(l-t)dt 
Jo 

= / (b"(uH + t(u h -u H ))(uh-u H ) 2 ,Xh)dt. 
Jo 

By Theorem 3.2, we have u h) u H G [u, u]. Therefore, we have the following estimate on 

R(u H ,Uh,Xh) ' 

\R{u H) u h) Xh)\ < C\\u H - u h \\l i2p \\xh\\o, q , (4.1) 

for any p,q>l with - + - = 1. 

Lemma 4.1. Let b satisfy the Assumption 2.2 and 3.4, and Th-> Th satisfy the Assump- 
tion 3.1. Let uh G Vh and Uh G Vh be the exact solutions to (3.1) on Th and Th 
respectively, and u h G Vh be the approximated solution obtained by Algorithm 1. Then, 
we have the following estimate 

UK -^111 % IK -uh\\Ia- (4-2) 
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Proof. From Algorithm 1, we have 

(F'(u H )(u h -u h ),Xh) = (F\u H )(u h -u H ),Xh) + (F(u H ),Xh) 

= -R{u H) u h) Xh), Vxh e 14. 

Then by taking Xh = u h — u h £ Vh in the above equality, we obtain that 

II ?Xfc — i^||| 2 = a(uh — u h , Uh — u h ) 

= (F'(u H )xh, Xh) ~ {b\u H )Xh, Xh) 

< -R(u H) u h) u h -u h ) 

< \\u h - u H \\l j2p \\u h - u h \\ 0jq , 

where in the third step we used the Assumption 3.4, and in the last step we used (4.1). In 
particular, if we pick p = q = 2, the conclusion then follows by Poincare inequality on 

11^-1^110,2. □ 

Lemma 4.1 suggests that we will need the L 4 error estimates \\u — u h \\o A . For this 
purpose, we make use of the following Ladyzhenskaya's inequalities: 

Lemma 4.2 ([19, Lemma 1-2]). For any v G Hq(Q), it holds 

||v||o,4<^||v|||, 2 ||Vv|| % d = 2- (4.3) 

and 

||v||o,4^V^||v||J 2 ||Vi;||| 2 , d = 3. (4.4) 

Recall that we assume the solution to the original nonlinear problem (2.1) satisfies 
u G Hq(Q) n H s (Qi U Q 2 ), an d the dual linear problem (3.9) has the regularity w G 
flo (fl) n i7 r (^i U Q 2 ) for some 5, r > 1. We let t = min{s, r} - 1 as defined in (3.17). 
As a corollary of Lemma 4.2, we obtain the L 4 error estimate: 

Corollary 4.3. Let b satisfy the Assumptions 2.2 and 3.4, and Th satisfy the Assump- 
tion 3.1. Let u G Hq(Q) n H s (Qi U Q 2 ) with s > 1 be the solution to (2.5), and Uh G Vh 
be the solution to (3.1). Suppose that the dual problem (3.9) satisfies the 7 -regularity 
(3.10) for some r > 1. Then the following error estimates hold: 

\\ u ~ u h\\oA < h^ 8 ~ 1) \\u\\ H s(Q 1 uQ 2 ), d = 2; (4.5) 

and 

\\ u ~ ^ II 0,4 % ^ i+(s_1) lkl|if s (^iu^ 2 ), d = 3, (4.6) 

where t = min{<s, r} — 1. 

Proof. The proof is simply a combination of Lemma 4.2 and the quasi-optimal error 
estimate (3.7) in Theorem 3.5 and the L 2 error estimate in (3.17). When d = 2, by (4.3) 

we have 

1 1 

\\u - u h \\ 0A ^ V2||ix - u h \\^ 2 \\^( u - ^)|| 2 52 < h t/2 \\V(u - u h )\\ 0i2 . 

Thus, when d = 2 we obtain 

Ik - ^/i 1 1 0,4 ^ ^^^"^IhllH-cniuna) 
Similarly, when d = 3 by (4.4) and (3.17) we obtain 

lh-^||o,4 < Ch*\\V(u - u h )\\ 0i2 . 
The conclusion then follows from the conclusion of Theorem 3.5. D 

Finally, we obtain the following main result. 



TWO-GRID METHODS FOR SEMILINEAR INTERFACE PROBLEMS 



13 



Theorem 4.4. Let b satisfy the Assumptions 2.2 and 3.4, and 7h, Th satisfy the Assump- 
tion 3.1. Let u e Hq(Q) fl H s (Qi U f^) with s > 1 be the solution to (2.5), and u h be 
the solution to Algorithm 1. Suppose that the dual problem (3.9) satisfies the r-regularity 
(3.10) for some r > 1. We have the following estimates 

h-u h \\\ < (h 3 - 1 + h^-V) 

and 



Mlif s (ftiuQ 2 )> 



d = 2, 



(4.7) 



\\u ~ u h \\\ < (/1 s - 1 + H-^^-^)\\u\\ Hs{niUn2h d = 3, 



(4.8) 



where t = min{s, r} — 1 > 0. 
Pro*?/ By triangle inequality, we have 



11^ — u III < \\\u — uA\ + | \\Uh — u 



< 

< 



\\u 



\\u 



Uh\\\ + \\u h - ^if 1 1 0,4 

^11 + \\ U - ^llo,4 + 11^ - ^if|lo,4- 



(4.9) 



The first term in the right hand side of (4.9) has been estimated in Theorem 3.5. Thus 
the conclusions immediately follow by applying Corollary 4.3 to the last two terms in 
(4.9). ^ ^ 

s-l 

Remark 4.5. Based on Theorem 4.4, we may choose H ^ h^ 2 ^- 1 ) in 2D and H < 

s-l 

j l t/2+2( s -i) ( n 3D b u t still achieve quasi-optimal error estimate. In particularly, if the 
linear dual problem (3.9) has the same or more regularity than the primal nonlinear 
problem (2.1), i.e., r ^ s > 1, then t = s — 1. So, in 2D case we may choose H ^ h 1 ^ 
and H ^ h 2 / 5 in 3D. 

5. Numerical Experiments 

In this section, we present some numerical experiments to justify the theories. Here 
we consider solving the following semilinear equation 

-V • (DVu) + u 11 = 10005 , u\ m = 0, 

where the diffusion coefficient D = 1000 inside [—1/2, 1/2] 2 and 1 outside, and S is the 
delta function on origin. See Figure 1 for the solution of this equation. Figure 2 and 3 




Figure 1 . Solution to the 2D semilinear interface problem. 
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Figure 2. Errors in i/^norm. 
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Figure 3. Errors in L 4 norm. 



show the H 1 and L 4 errors, respectively. Figure 4 shows the comparison of the exact H 1 
error with the error of the two-grid solution produced by the Algorithm 1 . Here, the mesh 
size H of the coarse grid problem is chosen to be closest to the theoretical ones obtained 
from Theorem 4.4 if not exactly the same. As we can see from this figure, the two-grid 
solution is very close to the exact solution. Therefore, by the appropriately choice of 
the coarse problem, solving the nonlinear problem could be reduced to solving a linear 
problem on the fine mesh without loss of accuracy. Note that the linearized problem on 
the fine mesh could be solved efficiently by multilevel preconditioning techniques, even 
in the presence of large jump coefficients (cf. [35]). In this way, we reduced greatly the 
overall computational cost for solving the nonlinear PDEs. 
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Figure 4. Two-grid Error. 



6. Conclusion and Extension 

In this article we considered a two-grid finite element method for solving semilinear 
interface problems in d space dimensions, for d = 2 or d = 3. We first described in some 
detail the target problem class with discontinuous diffusion coefficients, which included 
critical (and subcritical) nonlinearity examples, as well problems containing supercriti- 
cal nonlinearity (such as the Poisson-Boltzmann equation and the semi-conductor device 
modeling equations). We then developed a basic quasi-optimal a priori error estimate 
for Galerkin approximations. In the critical and subcritical cases, we follow [4] and 
control the nonlinearity using only pointwise control of the continuous solution and a 
local Lipschitz property, rather than through pointwise control of the discrete solution; 
this eliminates the requirement that the discrete solution satisfy a discrete form of the 
maximum principle, hence eliminating the need for restrictive angle conditions in the 
underlying mesh. However, the supercritical case continues to require such conditions 
in order to control the nonlinearity. We then designed a two-grid algorithm consisting 
of a coarse grid solver for the original nonlinear problem, and a fine grid solver for a 
linearized problem. We analyzed the quality of approximations generated by the algo- 
rithm, and proved that the coarse grid may be taken to have much larger elements than 
the fine grid, and yet one can still obtain approximation quality that is asymptotically as 
good as solving the original nonlinear problem on the fine mesh. The included numerical 
experiments support our theoretical results. 

The algorithm we described, and its analysis in this article, combined four sets of tools: 
the work of Xu and Zhou on two-grid algorithms for semilinear problems [32, 33, 34]; the 
recent results for linear interface problems due to Li, Melenk, Wohlmuth, and Zou [21]; 
recent work on the Poisson-Boltzmann equation [10, 14]; and recent results on a priori 
estimates for semilinear problems, including estimates without angle conditions in the 
case of sub- and super-critical nonlinearity [4]. Although the algorithm described in 
this paper is applicable to general coupled nonlinear elliptic systems, our reliance on 
tools developed for scalar linear and semilinear problems restricts the validity of the 
theoretical results to the class of semilinear problems described in §2. In future work we 
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will consider the case of coupled systems of scalar semilinear PDE from this class, as 
well as more general nonlinear elliptic systems. 

To simplify the presentation and keep the paper focused, we assumed that the triangu- 
lations resolve the interface. For general interface T, namely, T can not be resolved by 
the triangulation, we could use the concept of "5-resolved triangulation" (cf. [21]). The 
results in this article could be generalized in a direct way if the triangulation satisfies the 
"5-resolved". Without significant technical modifications to the results in the article, we 
could also relax the Local Monotonicity Assumption 3.4 to the following: 

6'(0 > -Ax, 

where Ai is the smallest eigenvalue of the operator —V • (DV-). 
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